# -*- coding: utf-8 -*-
"""Convert coordination in China
including GCJ-02, BD-09, WGS-84
@author: 红领巾1994
original code: https://www.cnblogs.com/giserjobs/p/12291291.html
@author: wandergis
original code: https://github.com/wandergis/coordTransform_py
Rewrite by SidChen"""


import math
import xlrd
from xlutils.copy import copy
import os


def file_name(file_dir):
    """Find goal file path"""
    l = []
    # root: current path, dirs
    for root, dirs, files in os.walk(file_dir):
        for file in files:
            if os.path.splitext(file)[1] in ('.xls', '.xlsx'):
                l.append(os.path.join(root, file))
    return l


def bd09_to_gcj02(bd_lon, bd_lat):
    x_pi = 3.14159265358979324 * 3000.0 / 180.0
    x = bd_lon - 0.0065
    y = bd_lat - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    # minus based on test
    return [gg_lng-0.012, gg_lat-0.012]

'''
def gcj02_to_bd09(lng, lat):
    x_PI = 3.14159265358979324 * 3000.0 / 180.0
    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_PI)
    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_PI)
    bd_lng = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    # minus based on test
    return [bd_lng-0.006264, bd_lat-0.001279]


def wgs84_to_gcj02(lng, lat):
    PI = 3.1415926535897932384626
    ee = 0.00669342162296594323
    a = 6378245.0
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * PI
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * PI)
    mglat = lat + dlat
    mglng = lng + dlng
    return [mglng, mglat]


def gcj02_to_wgs84(lng,lat):
    PI = 3.1415926535897932384626
    ee = 0.00669342162296594323
    a = 6378245.0
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * PI
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * PI)
    mglat = lat + dlat
    mglng = lng + dlng
    return [str(lng * 2 - mglng), str(lat * 2 - mglat)]


def transformlat(lng, lat):
    PI = 3.1415926535897932384626
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * \
          lat + 0.1 * lng * lat + 0.2 * math.sqrt(abs(lng))
    ret += (20.0 * math.sin(6.0 * lng * PI) + 20.0 *
            math.sin(2.0 * lng * PI)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * PI) + 40.0 *
            math.sin(lat / 3.0 * PI)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * PI) + 320 *
            math.sin(lat * PI / 30.0)) * 2.0 / 3.0
    return ret


def transformlng(lng, lat):
    PI = 3.1415926535897932384626
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(abs(lng))
    ret += (20.0 * math.sin(6.0 * lng * PI) + 20.0 *
            math.sin(2.0 * lng * PI)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * PI) + 40.0 *
            math.sin(lng / 3.0 * PI)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * PI) + 300.0 *
            math.sin(lng / 30.0 * PI)) * 2.0 / 3.0
    return ret


def bd09_to_wgs84(lng, lat):
    lng1, lat1 = bd09_to_gcj02(lng, lat)
    return gcj02_to_wgs84(lng1, lat1)


def wgs84_to_bd09(lng, lat):
    lng1, lat1 = wgs84_to_gcj02(lng, lat)
    return gcj02_to_bd09(lng1, lat1)'''

x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方


def gcj02_to_bd09(lng, lat):
    """
    火星坐标系(GCJ-02)转百度坐标系(BD-09)
    谷歌、高德——>百度
    :param lng:火星坐标经度
    :param lat:火星坐标纬度
    :return:
    """
    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)
    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)
    bd_lng = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    return [bd_lng, bd_lat]


def bd09_to_gcj02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    x = bd_lon - 0.0065
    y = bd_lat - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return [gg_lng, gg_lat]


def wgs84_to_gcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [mglng, mglat]


def gcj02_to_wgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [lng * 2 - mglng, lat * 2 - mglat]


def bd09_to_wgs84(bd_lon, bd_lat):
    lon, lat = bd09_to_gcj02(bd_lon, bd_lat)
    return gcj02_to_wgs84(lon, lat)


def wgs84_to_bd09(lon, lat):
    lon, lat = wgs84_to_gcj02(lon, lat)
    return gcj02_to_bd09(lon, lat)


def _transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
          0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def _transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)


def check_mode(mm):
    for k in range(0, 6):
        if not mm in ['11', '12', '21', '22', '31', '32']:
            mm = input('转换方式输入错误，请重新输入\n')
        else:
            break
    return mm


# end remark
def end():
    check = input_wait('\n* 检查脚本更新情况，稍候片刻…………\n（无需更新，可按Enter键退出）\n', 'checking......')
    if len(check) > 0:
        update_check(url_update, progress)
    print('\n'
          '\n****************************'
          '\n*坐标已全部转换，Enter键退出*'
          '\n****************************\n\n')
    input()
    exit()


def DownloadHtml(url):
    import string
    from urllib.parse import quote
    import requests
    import time
    from datetime import datetime

    time.sleep(5)
    # @Readme : 反爬之headers的伪装https://blog.csdn.net/weixin_44517301/article/details/103393145
    # 对于检测Headers的反爬虫
    from fake_useragent import UserAgent
    # 实例化，需要联网但是网站不太稳定-可能耗时会长一些
    ua = UserAgent()
    # 在工作中常用的则是ua.random方式
    headers = {
        # 伪装
        'User-Agent': ua.random
    }
    # 请求
    # 解决请求路径中含义中文或特殊字符
    url = quote(url, safe=string.printable)
    try:
        request = requests.get(url=url, headers=headers, timeout=(5, 27))
        html = request.text
        request.close()
        return html
    except:
        print('{}: 网页响应错误，30s后自动无限重试……\n（如果疯狂重试不成功，请手动关闭脚本，确认账号后再运行）'.format(datetime.now()))
        time.sleep(30)
        DownloadHtml(url)


def update_check(url, progress):
    import re
    from datetime import datetime
    import webbrowser
    from time import sleep

    # check update
    progress_ud = {}
    # progress_ud = {'name':[version_new, version_old, interval]}
    html = DownloadHtml(url)
    for i in progress:
        print(f'正在检查更新：{i}')
        # progress = {'name': ['version_old', 'download_link']}
        try:
            date = re.search(f'{i}.*"Modified.": "', html, re.S).span()[1]
            date_update = html[date:date+10].replace('-', '')
            version_new = datetime.strptime(date_update, '%Y%m%d')
            version_old = datetime.strptime(progress[i][0], '%Y%m%d')
            interval = version_new-version_old
            if interval.days > 0:
                progress_ud[f'{i}'] = (date_update, progress[i][0], interval)
                print(f'# 有更新:{date_update} ')
        except:
            print('* 未知错误')
    if len(progress_ud) > 0:
        load = input_wait('需要调起默认浏览器下载最新脚本文件吗？\n(默认需要，Enter退出)\n', 'Yes')
        if len(load) > 0:
            for file in progress_ud:
                # start the default browser to download the new version
                webbrowser.open(progress[file][1], new=2)
    else:
        print('* 无更新内容')


def input_wait(remind, defaut):
    """input or leap over"""
    # CODE from: https://blog.csdn.net/weixin_39858881/article/details/107152961
    import func_timeout

    @func_timeout.func_set_timeout(5)
    def askChoice():
        return input(remind)

    try:
        s = askChoice()
    except func_timeout.exceptions.FunctionTimedOut as e:
        s = defaut
    print('\n', s)
    return s


if __name__ == '__main__':
    print('脚本用于中国国内常用坐标系转换（WGS84，GCJ-02，BD-09）')
    # check and download the updated progress
    progress = {'Convert.exe': ['20210419',
                                'https://chdeducn-my.sharepoint.com/personal/201541020106_chd_edu_cn/_layouts/15/download.aspx?UniqueId=77e1a538%2D4284%2D4ca2%2Db9ef%2Dbcd62711c589'
                                '']
        , 'UpdateLogs.txt': ['20210419',
                             'https://chdeducn-my.sharepoint.com/personal/201541020106_chd_edu_cn/_layouts/15/download.aspx?UniqueId=480d6699%2Dde0b%2D476a%2Dba15%2Dfcc102c01ebf'
                             '']
                }
    # url must be the direct sharing link if it's from OneDrive where you can check latest update date
    url_update = 'https://chdeducn-my.sharepoint.com/:f:/g/personal/201541020106_chd_edu_cn/EpRLgdPCtklNgH-RfCoE_KkBZn0ZvgQHjAuqkX0iI1LLmQ?e=t3cc2t' \
                 ''
    for i in progress:
        print(f'* {i}，版本号{progress[i][0]}，')
    print('*** 有限技术支持：sidchen0 @ qq.com ***\n\n')

    # find current file path
    print('\n#############读取目标文件…………')
    p1 = os.getcwd()
    # find all excel path
    p = file_name(p1)
    pl = len(p)
    for i in range(0, pl):
        print('已找到\n{}:{}\n'.format(i, p[i]))
    # save work mode
    mo = {'11': 'bd09', '12': 'WGS84', '21': 'GCJ02', '22': 'WGS84', '31': 'GCJ02', '32': 'bd09'}
    print('\n#############选择转换模式')
    wm = []
    if pl == 0:
        print('当前目录及子文件夹下均未发现XLS/XLSX表格，任意键退出')
        input()
        exit()
    # one table
    elif pl == 1:
        m = input('请输入对应数字，选择坐标转换方式(高德，腾讯，谷歌中国为火星坐标系GCJ02)\n'
                  '11: 火星To百度\n12：火星ToWGS84\n21：百度To火星\n22：百度ToWGS84\n31：WGS84To火星\n32：WGS84To百度\n')
        m = check_mode(m)
        wm.append({'path': p[0], 'mode': m})
    # more than one table
    else:
        m2 = input('找到不止一个表格文件，请选择坐标转换模式\n1：相同转换模式\n2：不同转换模式\n')
        # one mode
        if m2 == '1':
            m = input('请输入对应数字，选择坐标转换方式(高德，腾讯，谷歌中国为火星坐标系GCJ02)\n'
                      '11: 火星To百度\n12：火星ToWGS84\n21：百度To火星\n22：百度ToWGS84\n31：WGS84To火星\n32：WGS84To百度\n')
            m = check_mode(m)
            for i in p:
                wm.append({'path': i, 'mode': m})
        # different modes
        elif m2 == '2':
            print('请输入表格编号及相应的转换模式，格式“表格编号,模式”（逐一确定或范围确定，英文逗号），如“1,11”或“1-10,22”\n'
                  '11: 火星To百度\n12：火星ToWGS84\n21：百度To火星\n22：百度ToWGS84\n31：WGS84To火星\n32：WGS84To百度\n'
                  'e，回车完成输入\n')
            for i in range(0, pl):
                p2 = input()
                try:
                    if p2 != 'e':
                        if not '-' in p2:
                            p2, m = p2.split(',')
                            m = check_mode(m)
                            wm.append({'path': p[int(p2)], 'mode': m})
                        else:
                            p2, m = p2.split(',')
                            m = check_mode(m)
                            p3, p4 = p2.split('-')
                            for j in range(int(p3), int(p4)):
                                wm.append({'path': p[int(j)], 'mode': m})
                    else:
                        break
                except:
                    print('!处理出错，确认是英文输入状态吗？\n（Enter退出，重启程序）')
                    input()
                    exit()
        else:
            print('输入错误，请任意键退出后重启脚本')
            input()
            exit()
    # read data
    # find coordination columns
    print('\n#############识别经纬度列')
    clgt = ['lgt', 'lng', 'lon', '经度']
    clat = ['lat', '纬度']
    print('以下列名称将被自动识别', clgt, clat)
    # complete work info
    wm2 = []
    for i in range(0, len(wm)):
        f = wm[i]['path']
        data = xlrd.open_workbook(f)
        table = data.sheets()[0]
        column_n = table.ncols
        # find coordination columns automatically
        h1 = -1
        h2 = -1
        for j in range(0, column_n):
            h = table.cell_value(0, j)
            if h in clgt:
                h1 = j
            elif h in clat:
                h2 = j
        if h1 >= 0 and h2 >= 0:
            print('已自动识别该表格的经纬度列：'+f)
        else:
            # find coordination columns physically
            print('无法找到给表格的经纬度列，请手动确认：'+f)
            for j in range(0, column_n):
                print(j, table.cell_value(0, j), table.cell_value(1, j))
            h = input('请输入经纬度列，格式“经度,纬度”(英文),如“1,4”\n')
            h1, h2 = h.split(',')
        wm2.append({'path': f, 'mode': wm[i]['mode'], 'lgt': h1, 'lat': h2})
    # convert and write in
    print('\n#############开始转换坐标')
    for i in range(0, len(wm2)):
        f2 = str(wm2[i]['path'])
        data2 = xlrd.open_workbook(f2)
        table2 = data2.sheets()[0]
        row_n = table2.nrows
        c = table2.ncols
        mode = wm[i]['mode']
        n1 = int(wm2[i]['lgt'])
        n2 = int(wm2[i]['lat'])
        # convert
        cd = []
        for j in range(0, row_n-1):
            lgt1 = float(table2.cell_value(j+1, n1))
            lat1 = float(table2.cell_value(j+1, n2))
            # convert
            if mode == 11:
                lgt2, lat2 = gcj02_to_bd09(lgt1, lat1)
            elif mode == 12:
                lgt2, lat2 = gcj02_to_wgs84(lgt1, lat1)
            elif mode == 21:
                lgt2, lat2 = bd09_to_gcj02(lgt1, lat1)
            elif mode == 22:
                lgt2, lat2 = bd09_to_wgs84(lgt1, lat1)
            elif mode == 31:
                lgt2, lat2 = wgs84_to_gcj02(lgt1, lat1)
            else:
                lgt2, lat2 = wgs84_to_bd09(lgt1, lat1)
            cd.append((lgt2, lat2))
        # write in
        cdata = copy(data2)
        ctable = cdata.get_sheet(0)
        ctable.write(0, c, 'lgt_'+mo[mode])
        ctable.write(0, c+1, 'lat_'+mo[mode])
        for j in range(0, row_n - 1):
            ctable.write(j + 1, c, cd[j][0])
            ctable.write(j + 1, c+1, cd[j][1])
        f3 = f2.rstrip('x')
        cdata.save(f3)
        print('\n坐标转换完成，已写入:' + f3)
    end()
